I 



SOME NUMERICAL COMPARISONS OF 
THREE-BODY TRAJECTORIES WITH 
PATCHED TWO-BODY TRAJECTORIES 
FOR LOW THRUST ROCKETS 

by William C. Strack 
Lewis Research Center 

t * 7 „ 

Cleveland \ Ohio ' :\ 

NATIONAL AERONAUTICS AND SPACE ADMINISTRATION % WASHINGTON, D. C. % MAY 1968 


TECH LIBRARY KAFB, NM 



tech library kafb, nm 


Q13151b 


SOME NUMERICAL COMPARISONS OF THREE-BODY TRAJECTORIES WITH 
PATCHED TWO- BODY TRAJECTORIES FOR LOW THRUST ROCKETS 

By William C . Strack 

Lewis Research Center 
Cleveland, Ohio 


NATIONAL AERONAUTICS AND SPACE ADMINISTRATION 


For sale by the Clearinghouse for Federal Scientific and Technical Information 
Springfield, Virginia 22151 — CFSTI price $3.00 




ABSTRACT 


The n-body problem for low-thrust vehicles with optimal thrust control is formulated 
for the case of constant thrust with coast. A successive two-body approximation to the 
n-body problem is also developed under the assumption of tangential thrust during the 
plane toe entric leg. This approximation involves one arbitrary parameter - the planeto- 
centric radius at the patch point - whose best value can be determined by comparison with 
n-body solutions. Numerical solutions of several three- body problems indicate that the 
patch radius should be about 300 Earth radii for either a circular or a parabolic initial 
Earth orbit. 



SOME NUMERICAL COMPARISONS OF THREE-BODY TRAJECTORIES WITH 
PATCHED TWO-BODY TRAJECTORIES FOR LOW THRUST ROCKETS 

by William C. Strack 
Lewis Research Center 

SUMMARY 

The n-body problem for low-thrust vehicles with optimal thrust control is formulated 
for the case of constant thrust with coast. Numerically calculated solutions are then pre- 
sented for several transfers between an initial low- Earth orbit (circular or parabolic) and 
the orbit of Mars. A successive two-body approximation to the n-body problem is devel- 
oped under the assumption of tangential thrust during the planetocentric leg. This approx- 
imation involves one arbitrary parameter - the planetocentric radius at the patch point - 
the best value of which can be determined by comparison with the n-body solutions. This 
value is found to be about 300 Earth radii for the circular initial Earth orbit case and for 
three different parabolic initial Earth orbit cases. If zero thrust is assumed for the 
planetocentric leg in the parabolic orbit case, the best value of the patch radius is between 
100 and 160 Earth radii. 


INTRODUCTION 

The prospect of using low acceleration propulsion devices in interplanetary vehicles 
has led to the development of elaborate mathematical techniques (e. g. , calculus of vari- 
ations schemes) that are used to optimize the thrust program along vehicle trajectories. 
These analyses are customarily carried out under the assumption of two-body motion and 
have proved quite useful in the solution of two broad classes of problems involving (1) two 
bodies only and (2) more than two bodies. The usual case involves a transfer from Earth 
to another planet. When the planetocentric legs of the trajectory are to be performed with 
high-thrust systems they are often ignored. The simplified problem then involves motion 
in a single gravitational field. The usual assumption is that low thrust commences in 
heliocentric space with a certain hyperbolic velocity (often zero) added to Earth's orbital 
velocity to determine the initial heliocentric velocity. The time required to escape from 



Earth orbit is usually neglected and the initial heliocentric radius vector is assumed to be 
that of the Earth. 

If the plane toe entric maneuvers are performed with low-thrust systems a patching 
technique is ordinarily used wherein only two-body motion is assumed at any given time. 
Quite often tangential thrust spirals to escape energy are calculated with approximate 
analytical expressions. If the vehicle’s planetoc entric patch radius and velocity s se as- 
sumed to be negligible in comparison to the Earth's heliocentric radius and velocity, then 
only the escape time and the propellant consumption during escape are needed to com- 
pletely determine the heliocentric conditions at the patch point. The same approach is 
often used when the rocket is assumed to be on a parabolic coast path while in planeto- 
c entric space. 

Several schemes following an approach commonly used in high- thrust trajectory 
studies have appeared recently (refs. 1 and 2) that introduce more realism into the models 
by relating the initial heliocentric conditions to both: (1) the planet's position and velocity 
and (2) the terminal planetoc entric radius and velocity of the rocket. The switch from 
planetoc entric to heliocentric coordinates takes place at a particular planetoc entric radius 
called the patch radius in this report. The initial heliocentric conditions (mass, time, 
radius, and velocity) are then dependent upon the patch radius, the planetoc entric thrust 
mode, and the size and shape of the initial planetoc entric orbit. The planetoc entric es- 
cape path can be determined either by precise numerical calculations or by approximate 
analytical expressions. In either case, the right patch radius and the proper spot on the 
patching sphere must be determined. The latter is easily taken care of by the transver- 
sality condition of the calculus of variations. The "right" patch radius is not so easily 
determined since the radius used in high-thrust studies is not necessarily applicable. If 
it is picked "properly, " the approximate two- body calculation will yield the same mass 
ratio as the n-body calculation. The corresponding trajectories, however, will always 
differ slightly. 

To date, no results have been presented comparing these approximate patching 
schemes with n-body variational solutions because the n-body solutions are difficult to 
obtain and require a large expenditure of computer time. The purpose of this report is 
to (1) present several three-body solutions, and (2) compare the approximate two-body 
patching technique results with these more exact solutions. 

Three-body solutions are generated for low thrust transfers from low Earth orbits to 
the mean heliocentric orbit of Mars. These solutions are for a simple problem model 
involving three bodies: the rocket, Sun, and Earth. This simplification in the sample 
problems allows attention to be focused on the main issue - trajectories that pass from 
one dominant gravitational field into another. The thrust vector is optimally directed and 
of constant magnitude. Coast phases are permitted. Two cases are considered: (1) ini- 
tial circular orbit about Earth and (2) initial parabolic orbit about Earth. These cases 
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correspond to escapes utilizing (1) all low-thrust and (2) high thrust followed by low 
thrust. The corresponding patched two- body calculations are carried out under the same 
rules except that (1) tangential thrust is assumed for the Earth-centered portion of the 
trajectory and (2) the two-body assumption is imposed. The tangential thrust assumption 
is not critical since tangential thrust is very nearly optimum for low-thrust Earth escapes 
(ref. 3). Some solutions are also presented for the parabolic orbit case under the as- 
sumption that the vehicle coasts to the patching point. 


ANALYSIS 

The problem posed is: Given the initial state of a rocket in Earth orbit and a helio- 
centric final state, determine the low-thrust transfer trajectory that minimizes the pro- 
pellant consumption. The accurate solution of such a problem is difficult and quite time 
consuming, and no results have appeared in the literature. A second question, therefore, 
arises: How can the actual solutions be closely approximated? This question will be dis- 
cussed later. 


N-Body Solutions 

Statemen t of problem. - The rocket is assumed to operate with constant exhaust ve- 
locity engines that may be freely turned on or off. Furthermore, the thrust magnitude 
remains constant when the engines are operating. The vector equation of motion of such 
a rocket in a central force field perturbed by n-2 other bodies is: 


n-1 


R = 
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(1) 


Figure 1 displays the geometry of the problem and all symbols are defined in appendix A. 
The center of the dominant gravitational body defines the origin of the coordinate system 
and the summation term vanishes for two-body motion. 

The necessary conditions for minimizing the propellant consumption are given by the 
Euler- Lagrange equations (ref. 4): 
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(2) 



where X is the Lagrange multiplier vector (also called the adjoint variable or primer 
vector). Rj, r^ and pj in equation (2) are identical to R, r, and p in equation (1). 
The choice of on or off mode of engine operation is determined by: 


m = 0 if j^c | X | - am < oj 

(5) 

™ = “max if [ c M am > oj 

(6) 


Equations (1) to (6) comprise a set of second order nonlinear differential equations 
that must be numerically integrated to yield solutions. Equations (1) and (2) are equiv- 
alent to the following first order set of differential equations: 



Just as in the case of two-body motion, the transversality condition (ref. 4) is: 
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-drn^ + |\ • dV + A • dR + adm - C dt]* = 0 


( 11 ) 


where o and f denote values taken at the beginning and end of the optimal control phase, 
respectively. This equation imposes boundary conditions in addition to those explicitly 
required in the problem definition. On a solution trajectory a differential vanishes at 
either end point in equation (11) if its corresponding variable has a specified value at that 
terminus, otherwise the differential is nonzero and its coefficient must vanish. In par- 
ticular, since m^ is not specified, equation (11) shows that Oj = 1. C is a constant of 
the motion and is given by 


C=A-V+A*R + crm = constant (12) 

The following additional vector constant of the motion is valid for the two-body problem 
only (ref. 5): 


A XR + RXA = a vector constant (13) 

Two types of discontinuities may occur along a trajectory. First, m may switch from 
m max to zero or vice-versa. The Weierstrass- Erdmann corner condition (ref. 4) stip- 
ulates that A, A, and a are continuous across such corners and it is easy to show that 
C is also continuous at such points. Second, in the case of n-bodies, a coordinate sys- 
tem change may take place at a certain planetocentric radius. Thus, R and V are sub- 
ject to discontinuities. At the translation point C is also discontinuous. 

Boundar y val ue problem. - The boundary conditions at the initial and final points are 
satisfied by successive iterations of the initial values of A and A. In practice, the ini- 
tial value of a may be arbitrarily chosen since the set of equations (4) to (10) is homo- 
geneous and the = 1 condition may be satisfied simply by rescaling the multipliers 
after the boundary- value problem is solved. Many iterative techniques may be used to 
satisfy the terminal boundary conditions. The technique employed in this study is a hy- 
brid combination consisting of (1) successive univariate searches when the error function 
is relatively large and (2) a multivariate Newton- Raphson scheme when the error function 
is relatively small. The error function is simply a weighted summation of the squared 
differences between the final desired and current boundary conditions. Rather than com- 
puting the required partial derivatives of the Newton- Raphson scheme by finite differ- 
encing, the partial derivatives were generated by integrating an additional set of equations 
simultaneously with each trajectory calculation. The advantages of this technique are that 
convergence of the Newton-Raphson scheme is facilitated since the resultant partial deri- 
vatives are more accurate, decisions about the increments in the initial values of A 
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and X used for the perturbation trajectories are avoided, and less computer time is re- 
quired for a solution. For example, inclusion of a third body in the problem makes so- 
lutions very difficult to obtain by a finite difference scheme but relatively easy with ana- 
lytical partial derivative generation. These additional differential equations come directly 
from the equations of motion and the Euler- Lagrange equations. Differentiating equa- 
tions (4) to (10) with respect to an arbitrary variable y i and inverting the order of differ- 
entiation, after simplifying, yields: 
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when m is constant, 9m/9 y. = 0 and equation (19) reduces to the trivial identity 0 = 0. 
With the assumption of constant thrust, m is constant everywhere except the points at 
which c | X | - am = 0 in which case 3m/9y^ is unbounded causing jump discontinuities in 
all equations containing a 9m/9 y^ term. The appendix of reference 6 shows that the 
magnitude of these jumps is given in general by: 


A 
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(a continuous function everywhere). Specifically, 
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( 22 ) 


(23) 


(24) 


(25) 


(26) 


It should be noted that these partial derivatives are continuous under an origin trans- 
lation. Also, the partial derivative equations (14) to (19) are insufficient by themselves 
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for the case of unspecified terminal time. Only fixed terminal time cases are considered 
here. When coast phases are permitted the y i are taken to be the initial values of X 
and X. When the problem does not allow coat phases, a is superfluous and the initial 
value of any one of the components of X or X may be arbitrarily selected. For this 
case, one of the y ^ is taken to be m since for fixed time problems m is unique but 
unknown. The initial values of the partial derivatives are functions of the initial boundary 
conditions. For example, if the values of V and R Q are fixed then all partial deriva- 
tives are initially zero except [9X/3xj 0 and [9X/ax] 0 which are unit diagonal matrices. 
But if V Q and R q are some function of x , then the initial values of 3V/3 y^ and 
3R/07^ are not all zero. Adding a fixed amount of velocity via high thrust propulsion to 
some reference velocity (e. g. , the Earth's orbital velocity) in an optimal direction is an 
example in point. This will be made more explicit later. 

The Newton- Raphs on iteration scheme can be represented mathematically by: 

’'new = ’'old + k < Ar > < 27 ) 

Ay = A“ 1 (AB) (28) 

where k is a damping factor whose value lies between 0 and 1, AB is the vector of re- 
sidual errors in the boundary conditions, and [a] is the square matrix of partial deriva- 
tives dB^/dy-. The elements of [a] are composed of a subset of the definite integrals of 
equations (14) to (19) or combinations of these integrals depending upon the form of the 
end conditions. If the end conditions are given by (point to point problem): 


> 

II 

(29) 

R f = R 

(30) 


where a bar denotes a desired value, then [a] consists simply of the elements (dV/dy^ 
and (SR/Sy^. For the flyby problem with fixed terminal distance but free velocity and 
path angle, the end conditions are: 


r f - r 

(31) 

o 

ii 

(32) 
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and the elements of [a] consist of (Sr/Sy^ and (3X/3y^. Here (3X/3 y^) is determined 
by integrating equation (16) and (3r/3y^ is determined by integrating equation (15) via 
the relation 


3 Y t r \ 3 yj 


(33) 


Obviously many other sets of end conditions are possible and each of these results in a 
different form for the elements of [a]. 

Most of the results presented later are stated in terms of the propellant expenditure. 
In this connection, the final mass value may be linearly corrected for residual errors in 
the terminal boundary conditions. In the simplest case the terminal boundary conditions 
are fixed values of the state variables V and R. Since X and X may be interpreted 
(ref. 7) as the partial derivatives of the negative of the final mass with respect to their 
corresponding problem variables, 


m^ s m^ + X • (V - V) + X 



(34) 


X f and Xj are assumed to be normalized by cr f = 1 as required by the transversality 
condition. 

Procedure. - The specific problem for which numerical results are later presented 
is the following: At time zero the Earth is on the X-axis of a heliocentric coordinate sys- 
tem and a space vehicle is at the perigee of an orbit about Earth. The perigee altitude is 
185 km and the vehicle’s orbit about Earth is either circular or parabolic. At the speci- 
fied final time the vehicle arrives at the orbit of Mars after traversing 225° (arbitrarily 
chosen) of total heliocentric angle. For simplicity, the orbits of Mars and Earth are as- 
sumed to be circular and coplanar with the vehicle's initial orbit about Earth. Table I 
contains a list of assumed values of problem constants and table II displays the combina- 
tions of trip time and initial thrust-weight ratio F/W picked for sample cases. 

Sketch (a) shows the geometry. The angular location of the initial point with respect 
to the X-axis 0 Q is free to be optimized. The optimization can be accomplished by 
satisfying the transversality condition (eq. (11)) or, alternatively, by employing any sort 
of simple search scheme. Specifically, the transversality condition requires that 

^XjV cos 0 +XgV sin 0 +X^r sin 0 - Xgr cos 0^=0 (35) 

where X^ and X 2 are the X and Y components of X and similarly for X^, Xg, and X. 
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The numerical calculations were carried out with the present version of the Lewis 
N-Body Code which uses a fourth order Runge-Kutta integration scheme with variable step 
size control (ref. 8). Test cases showed that solution values were quite insensitive to 
the plane toe entric radius at which the origin was translated from the Earth to the Sun. 
Over the range of 25 to 10 000 Earth radii, the final mass change was less than two units 
in the sixth significant figures. The particular radius used for data generation was 
145 Earth radii. Obtaining solutions for n-body problems is substantially more difficult 
than for two-body problems. The computer time alone is one to two orders of magnitude 
greater for n-body solutions compared to two- body solutions. The solution method used 
for this study is briefly summarized as follows: First, obtain the solution of a relatively 
simple problem that only approximates the actual problem. The simple problem is ob- 
tained by setting p 2 ~ ® so body effects are ignored and by using a tangential 

thrust program until r = r* before switching to optimal thrusting. Second, gradually 
increase to its real value by solving a series of problems with pg as a parameter. 
Third, reduce r* to the initial radius r^ n through another series of problems. Finally, 
when r = r^ n the solution obtained is the actual n-body solution of the original problem. 
When one such solution has been obtained, the precise solution of other problems can be 
found by solving a series of problems starting in the neighborhood of the solved problem 
and proceeding toward the desired one. The only real difficulty with the above procedure 
is that as r* approaches r^ n the trajectories become increasingly more sensitive to 
changes in the initial values of A and X. This situation becomes so severe that conver- 
gence is impossible for r* less than about 20 Earth radii due to the limited number of 
significant figures carried by the computer (8 in this study). Although distressing, this 
situation is not so bad as to prevent meaningful answers from being obtained. This is 
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true because, as was stated in the INTRODUCTION, tangential thrust is very nearly opti- 
mum for small r, and therefore, a plot of against r* can be extrapolated to r^ n 
quite accurately (i. e. , to 4 significant figures). 


Approximate Patched Two-Body Solutions 

Patched two-body approximations to the low- thrust n-body problem are popular 
because of their computational simplicity. The assumption usually made is that the 
rocket coasts or is tangentially propelled to the patching point ignoring the Sun's gravi- 
tational effect. From that point on the rocket motion is calculated under the influence of 
no other gravitational body except the Sun. The equations of motion are much simpler 
and can be solved analytically for coast phases. When low-thrust propulsion commences 
from a circular planetocentric orbit, tangential thrust is most often assumed because it 
is very nearly optimal (ref. 3) up to the escape condition. When the initial planetocentric 
orbit is parabolic or hyperbolic, coast flight is usually assumed for the planetocentric leg, 
although tangential thrust is also sometimes assumed. Furthermore, the tangentially 
powered motion in planetocentric space is almost always represented with one of many 
analytic methods (e. g. , refs. 11 and 12). Optimal thrust programming begins in helio- 
centric space. 

In this study tangential thrust is assumed for the planetocentric leg and optimal thrust 
control is used for the heliocentric leg. Furthermore, the tangential thrust phase is 
numerically integrated rather than approximated in order to avoid results that are de- 
pendent upon specific approximations. In some cases the time optimal control policy 
calls for a coast phase at the end of the planetocentric leg. This situation can be detected 
by the presence of an initial coast phase on the heliocentric leg and in this case an itera- 
tion on the planetocentric coast duration is required to determine the best two-body 
patched trajectory. 

As in the n-body problem, the initial heliocentric position and velocity of the vehicle 
are obtained by vectorially adding the vehicle's planetocentric position and velocity at the 
patching point to the heliocentric position and velocity of the departure planet. However, 
since the initial planetocentric point may be anywhere on a sphere of radius r , the 
patching point on this sphere is also arbitrary. This arbitrariness is removed by re- 
quiring that the patch point be the optimal one. This in turn determines the optimal 
planetocentric starting point. Suppose that in general the position of the departure planet 
is Rp relative to a Cartesian coordinate system with the origin located at the center of 
the Sun (see sketch (b)). Let the planet's velocity be Vp and let the vehicle position and 
velocity at the patch point relative to a planet centered coordinate system (alined with the 
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Sun centered system) be R and V . Then, relative to the Sun: 

s s 


R o - R p + R s 


(36) 


V = V + V 
o v p v s 


(37) 


And from sketch (b), for the two-dimensional case, the components of R and V are 

s s 


x' - r cos 6„ 
s s 


(38) 


y' = r sin 0 
J s s 


(39) 


and 


x' - -v g sin (6 s - P 8 ) (40) 

y' = v g cos (0 g - /3 g ) (41) 

The transversality condition (eq. (11)) requires that 

V dV o + V dR o = ° < 42 ) 

Evaluating this expression with the aid of equations (36) to (41) and noting that R Vp, 

r , v and fi are all fixed, and solving for the optimum value of 6 yields 
s s s s 


I 
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(43) 


so, 


tan 0 s 


v g (-X 1 cos J3 g + Xg sin /3 g )+X 2 r s ^ ^ 
Vg^ sin /3 g +X 2 cos /3 g ) + E 



sin 0 g 


D 

i^D 2 + E 2 


(44) 


Convincing arguments that the plus sign is correct in equation (44) are given in refer- 
ence 2. Similar relations for the three dimensional case exist but are quite lengthy (two 
additional degrees of freedom are involved). As stated earlier, problems that have an 
initial optimal control point that is not fixed but lies on some curve or surface cause some 
of the analytical partial derivatives to have initial values dependent on the boundary sur- 
face function. The heliocentric leg of the two-body patched scheme outlined above is a 
case in point. From equations (36) and (37), 

\ 3R S 

= — - (45) 

L 3r * 




(46) 


These equations are solved for the cases y, = X, and y. = X, in appendix B. 

X 1 } U X X y U 

This completes the analysis of the n-body and patched two-body techniques used in 
this report. The remaining portion is devoted to several specific examples that compare 
three-body solutions with patched two-body solutions. 


RESULTS AND DISCUSSION 
Three-Body Solutions 

The sample problems all involve a transfer from an Earth orbit condition to the 
orbit of Mars. The main concern here is the computation of optimal thrust trajectories 
that essentially leave one gravitational field and enter another. All unnecessary com- 
plicating aspects of a real problem are avoided. Thus, the n-body problem is reduced to 
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a three-body problem by ignoring the minor perturbative effects of bodies other than the 
Sun and the Earth. Furthermore, the orbits of Earth and Mars are taken to be circular 
and coplanar with the vehicle's orbit. Three dimensional ephemeris data and additional 
bodies can just as well be included (i. e. , the computer code used contains such realism) 
but would serve no useful purpose here. 

The trajectories for case 1 and 2 of table II are diagrammed in figures 2 to 4. These 
cases are for 275-day/225° transfers with F/W = 10”^. Figure 2 shows the planeto- 
centric phase for the case of an initial parabolic orbit about Earth. As explained in the 
ANALYSIS section, a short tangential thrust phase is necessary at the start of the tra- 
jectory for computational convenience. The switch from tangential to optimal thrust con- 
trol takes place at 23 Earth radii. In Earth space the optimal angle- of -attack (angle from 
the velocity vector to the thrust vector) increases from 0. 5° at the initial point to only 
5. 2° at the origin translation point. Thus, it is quite unlikely that replacing the short 
tangential thrust phase with optimal control would cause any significant change in either 
the trajectory or the optimal angle -of -attack schedule. 

Figure 3 shows the planetocentric phase for an initial circular orbit about Earth. 

Only the last two turns of the spiral are shown and the optimal control policy begins at 
20 Earth radii. The angle -of -attack history is again nearly tangential but oscillates, 
starting at +4°, dipping to -13°, and then rising to -2. 5° at the origin translation point. 
This behavior has been noted before in the optimal Earth escape problem (ref. 3). It con- 
sists of minute oscillations of the thrust vector about the velocity vector until the last few 
turns of the spiral where the amplitude of the oscillation grows quite rapidly. 

The heliocentric portions of cases 1 and 2 are displayed in figure 4. Both trajectories 
are rather typical, having a mixture of power and coast phases. The parabolic orbit case 
has two coast phases but the circular orbit case has only one. The entire angle -of -attack 
histories are shown in figure 5 for both cases. The initial tangential thrust portion for 
the parabolic orbit case is so short that it is not detectable in the scale of this figure. 

Quite the reverse is true for the circular orbit case. Here the tangential thrust portion 
is seen to last about 67 days and represents the average angle- of -attack history of the 
true oscillatory solution. The optimal thrust angle history in inertial space is, of course, 
continuous across the origin translation. However, the discontinuity in the velocity vector 
at this point causes a discontinuity in the angle -of -attack history as seen in this figure. 

Mentioned earlier was the fact that a one-parameter sequence of solutions is obtained 
wherein the parameter is the radius of the switch point from tangential to optimal thrust 
control. Ideally, this sequence is terminated when the switching radius becomes equal to 
the initial radius thus eliminating the tangential thrust phase. Practically, numerical dif- 
ficulties prevent solutions from being obtained when the switching radius falls below about 
20 Earth radii. Nonetheless, such a terminated sequence allows the ideal solution to be 
accurately determined to more than four significant figures in the propellant expenditure 
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by extrapolation. Having determined the three-body solution (see table II for propellant 
fractions), the percent increase in propellant expenditure is plotted against the switching 
radius r* in figure 6. The two curves in figure 6 have discontinuous slopes at the origin 
translation radius because the tangential thrust vector is discontinuous there (i. e. , tan- 
gential thrust with respect to the Sun is different than tangential thrust with respect to 
Earth). This effect is much more pronounced in the parabolic orbit case because the ini- 
tial optimal control in heliocentric space is much farther from tangential than in the cir- 
cular orbit case (see fig. 5). In either case figure 6 suggests that if the rocket is in 
heliocentric space and not too far from the planet then tangential thrust with respect to 
the planetocentric velocity vector is better than tangential thrust with respect to the helio- 
centric velocity vector. The most important point though is that tangential thrusting in 
planetocentric space is not far from optimal in terms of propellant expenditure. That is, 
instead of insisting that r* = r^ n , r* may be as large as 145 Earth radii with less than a 
0. 12 percent penalty for the circular orbit case and less than 0. 01 percent penalty for the 
parabolic orbit case. The much smaller penalty for the parabolic orbit case is a result 
of the much smaller time period associated with these nonoptimal phases (e. g. , 1/2 day 
against 67 days for r* = 25 Earth radii). 

Cases 1 and 2 illustrate the essential features of the three-body solutions. Cases 3 
and 4 are presented mainly to show what effect changes in F/W and mission AV have 
on the comparison between three-body and patched two-body solutions. These compari- 
sons are given in the next section. 


Patched Two-Body Solutions 

Most of the patched two-body solutions used to approximate the three-body solutions 
consist of a tangentially propelled planetocentric leg followed by a heliocentric leg with 
optimal thrust control. Some results are also given for planetocentric coasting (instead 
of tangential thrusting) for the parabolic orbit case. In all cases the decision must be 
made of when to terminate the planetocentric leg and commence the heliocentric leg. The 
patch point is chosen to occur at a specified planetocentric radius r . This patch radius 
may, in fact, be so chosen that the correct propellant expenditure is calculated (i. e. , the 
propellant expenditure agrees with the three-body solution). Of course, the two-body and 
three-body trajectories will differ slightly, but that is not critical for most mission 
analysis studies. 

The fact that a correct patch radius exists is easy to see by examining the two 
limiting cases r g = r in and r g = °°. If r g = r^ n the Earth's gravitational field is lost 
and, therefore, (1) the propellant needed to escape the Earth is saved and (2) the propel- 
lant needed for the heliocentric leg is reduced. The heliocentric propellant is reduced 
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for two reasons. First, the initial heliocentric velocity will be most favorable since V 

s 

increases as r decreases. Second, the heliocentric travel angle and time will be a 

maximum thus allowing the most efficient trajectory. These effects taken together cause 

the required AV to be minimal. Thus, the two-body propellant mass is less than the 

three-body propellant mass if r is sufficiently small. * On the other hand, as r — °° 

s s 

the propellant mass increases until it approaches the initial mass or the plane toe entric 
time exceeds the total mission time. Clearly then, the two-body propellant mass exceeds 
the three-body propellant mass if r is sufficiently large. 

Tangential t hrust for plane toe entric leg. - The ratio of the two-body propellant mass 
to the three-body propellant mass is presented in figure 7 for cases 1 to 4 with tangential 
thrusting to r . The tangential portion is calculated exactly to determine the correct 

O 

time, velocity, and path angle at r . Equations (36) to (44) are then employed to deter- 

mine the optimum initial heliocentric position and velocity. Finally, equations (4) to (10) 

are integrated to determine the optimal heliocentric trajectory. The ratio of propellant 

masses increases as r increases in accordance with the reasons given in the previous 

paragraph. Figure 7 shows that at r g = 300 Earth radii the ratio is nearly unity for every 

case. This indicates that r « 300 Earth radii is the correct value for the two-body patch 

s 

radius. For r < 300 Earth radii not enough propellant is calculated by the two-body ap- 
s 

proximation and for r > 300 Earth radii too much propellant is calculated. Note also 

b 

that for r > 150 Earth radii the curves are fairly flat and the maximum propellant error 
in the two-body scheme for 150 < r < 600 Earth radii is about 3 percent. Larger er- 

b 

rors prevail at r g < 150 Earth radii. 

The case of an initial circular orbit about the Earth is the least sensitive of all cases 

investigated, r values between 250 and 600 Earth radii give rise to less than 1/2 per- 
s 

cent error in this particular case. Sometimes the patching radius for two-body schemes 
involving an escape spiral is taken to be that radius at which escape energy is attained. 

This radius is problem dependent and varies considerably with acceleration level. For 
the case considered here, escape occurs at 80 Earth radii and this gives rise to a -3 per- 
cent propellant error as noted on figure 7. The other three cases are all for an initial 
parabolic Earth orbit. Reducing the F/W from 10 to 0. 56x10 reduces the correct 

r from 310 to 285 Earth radii and increases the sensitivity of the error to r by about 
s s 

30 percent. Reducing the trip time from 275 days to 240 days (this increases the mission 

AV requirement by a factor of 2^) causes no change in the correct value of r g and only 

slightly decreases the sensitivity of the error to r . 

s 

Coastfor plane toe entric leg. - If the initial planetoc entric orbit is parabolic no thrust 

*In the extreme cases of very small mission AV and/or high F/W, the heliocentric 
powered phases can disappear entirely as r is reduced. In such a case the minimum 

b 

propellant case occurs at r g > r^ n> 
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is needed to reach the patching point. The vehicle may be allowed just to coast to this 
point. Assuming coast flight for the plane toe entric leg instead of tangential thrust for the 
parabolic initial orbit case leads to the results given on figure 8. The dot-dash curve is 
for the additional assumption V g = 0 and illustrates that this assumption causes the two- 
body approximation to be over 10 percent conservative for this particular problem. The 
other two curves show that the coast assumption alone increases the sensitivity of the er- 
ror to r and also reduces the correct value of r from about 300 Earth radii to 100 to 
s s 

160 Earth radii depending upon F/W. This marked increase in sensitivity and range of 
"correct r " is not surprising in view of the significant coast durations involved (e. g. , 

b 

23 days for r = 300 Earth radii). The tangential thrust assumption is definitely prefer- 
able to the coast assumption with respect to the sensitivity and range considerations. The 
trajectories will also be more accurate under the tangential thrust assumption since this 
mode of operation more closely approximates the three-body solution. 

Comparison of r s to analytic sphere of influence radius. - The classical sphere of 
influence radius is obtained by setting the ratio of the perturbative force to primary 
force in planetocentric coordinates equal to the same ratio in heliocentric coordinates 
(ref. 9). This gives rise to a value of 145 Earth radii for the Earth-Sun system. This 

value is roughly in agreement with the planetocentric coast results (100 < r < 160 Earth 

s 

radii) found in this study. But is is only 1/2 the value found under the tangential thrust 
assumption. 

Another sphere of influence definition (ref. 10) sets the Sun's perturbative force equal 
to the planetocentric force, computed in planetocentric coordinates. The formula (derived 
in ref. 10) for r using this definition is: 

b 


r = (planet's orbit radius)(i ? lanet mass \ 1/3 (47) 

\2 Sun mass / 

This leads to a sphere of influence radius of 270 Earth radii for the Earth-Sun system. 
This value compares quite well with the empirical values (280 < r g < 320 Earth radii) 
found here under the tangential thrust assumption. Extending the usefulness of equa- 
tion (47) to other planet-Sun systems is tempting because of the good correlation in the 
Earth-Sun case, but such an extension must be regarded as quite speculative at the pres- 
ent time. 

The natural suggestion based on the limited number of cases investigated is to em- 
ploy tangential thrust for the Earth-centered leg (since it agrees with the three-body 
power mode) and to use 300 Earth radii for the two- body patching radius r . 

b 
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CONCLUDING REMARKS 


The three-body variational solutions obtained provide reference data useful in the 
evaluation of approximation schemes. The optimal thrust direction in plane toe entric 
space is nearly tangential and this leads to the assumption of tangential thrust during the 
plane toe entric leg of the patched two-body scheme. Under this assumption, a value of 
300 Earth radii for the patch radius never causes more than 1/2 percent error in propel- 
lant consumption for each of the four cases. Of course, these four cases do not cover a 
wide enough range of problems to conclude that 300 Earth radii is a good value for most 
problems. Nonetheless, the fact that doubling the F/W, doubling the mission AV, or 
changing the initial Earth orbit from circular to parabolic causes such small disturbances 
in the error at this value leads one to believe that it may hold over a much greater range 
of problems than considered here. Other destinations, for instance, would probably not 
affect the Earth-centered leg significantly so that this type of problem alternation is 
unlikely to affect the answers found here. If a nontangential thrust program were speci- 
fied because of a system constraint (such as for solar- elec trie propulsion system), the 
results found here again would seem to be applicable. 

The plane toe entric coast assumption produces less uniformity in the correct value of 

r and significantly increases the propellant error at other r It also decreases the 
s s 

correct value of r to 100 to 160 Earth radii. Only parabolic initial Earth orbit cases 

b 

are examined here although the circular orbit case can also be done if the coast phase 
begins after escape energy is attained. In either case, a planetoc entric coast phase as- 
sumption simplifies the problem (permits closed form solutions to the planetoc entric tra- 
jectory) at the expense of accuracy. It should be noted here that if the three-body trajec- 
tories have a relatively large planetoc entric coast arc, instead of a thrusting arc, the 

right value of r may well be about 300 Earth radii instead of 100 to 160 Earth radii, 
s 

This conjecture is unresolved at the present time. 

It seems that the most accurate method of approximating the three-body solution is 
to numerically calculate both a powered planetoc entric leg and a heliocentric leg. It is 
quite possible that the need to numerically calculate the planetoc entric leg can be removed 
through the use of analytic approximations such as are found in references 11 and 12. 

This would certainly reduce the required computer time necessary to obtain solutions but 
still retain a high degree of accuracy. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, January 24, 1968, 

789-30-01-01-22. 
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APPENDIX A 


SYMBOLS 


A 

partial derivative matrix 

AV 

characteristic velocity incre- 

B 

residual error vector 


ment, m/s 

C 

first integral of Euler- Lagrange 
equation 

V P 

velocity vector of a planet rel- 
ative to the Sun, m 

c 

D 

rocket exhaust velocity, m/s 
defined by equation (43) 

V s 
v, v 

velocity vector of the rocket at 
the patching point, m/s 

magnitude of V and V , re- 

E 

defined by equation (43) 

D 

s 

spectively, m/s 

F/W 

initial thrust-weight ratio 

X, Y 

coordinate axes, origin at Sun 

g 

switching function defined by 


center 


equation (21) 

X',Y’ 

coordinate axes, origin at 

k 

damping factor 


planet center 

m 

mass of rocket, kg 

x’,y’ 

X’,Y' components of R m 

b 

R 

position vector of rocket, m 

0 

path angle 

R j 

■fVi 

vector from j gravitational 

body to rocket, m 

r 

vector of unknown initial con- 
dition variables, usually 

R P 

position vector of a planet 


composed of X Q and X Q 

relative to the Sun, m 

e 

central angle of rocket relative 

R s 

planetocentric position vector 


to the X-axis 


at which the two- body 
patching scheme is effected, 

X 

Lagrange multiplier vector as- 
sociated with velocity 


m 

X 

Lagrange multiplier vector as- 


radius at which optimal thrust 


sociated with position 


control begins, m 

W 

magnitude of vector X (X is 

r ’ r j’ r s 

magnitude of R, R., and R , 
J S 

respectively, m 


any vector) 

^3 

ii. 

gravitational constant of j 

S 3 

V, cr, or m 

3 2 

gravitational body, m /s 

T 

V 

unit vector in thrust direction 
rocket velocity vector, m/s 

a 

Lagrange multiplier associated 
with mass 
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Subscripts: 

f end of optimal control 

j-i_ 

i the i u component of a vector 

in start of complete trajectory 

j the j u gravitational body 

max maximum 


o start of optimal control 

s at the patching point 

Superscripts: 

desired value 

differentiation with respect to 
time 
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APPENDIX B 


INITIAL VALUES OF THE ANALYTICAL PARTIAL DERIVATIVES 

Derived here are the initial values of the analytical partial derivatives for the two- 
dimensional heliocentric leg of the two-body patched scheme. From equations (38) to 
(41), (45) and (46) the required derivatives have the initial values: 





(48) 


(49) 



cos (0 g 




sin (0 g 



where, from equation (43), 



(50) 


(51) 


(52) 


Specifically, if y is composed of the components of X and X , then from equation (43): 


8D 


dr 


l 



cos /3 g 


(53) 
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9E 


(54) 


37l 



= v. 


sin ^ 



5E 

dy 2 



COS 


(55) 


(56) 



Substituting equations (53) to (60) into equation (52): 



V2 +v s r s ( *l COS / 3 S +X 2 sin * B > 


D 2 + E 2 



Vl + v s r s (X l sin ^s"*2 cos *s> 


D 2 + E 2 


— ' o 


(57) 


(58) 


(59) 


(60) 


(61) 


(62) 
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(63) 



Vs ( - A l C0S ^s +X 2 sin ^s )+r s*2 


D 2 + E 2 



r s v s (X l sin + x 2 cos e*) + Vl 


D 2 + E 2 


(64) 


Finally, substituting equations (61) to (64) into each of the equations (48) to (51) leads to 
the sixteen equations that determine the initial values of the analytical partial derivatives 
3R /dr i and dV/dy,. These equations are not listed here in their completely written out 
form because of their length and number. 

The other partial derivatives have the initial value: 


dx i ax t 
ax. ax. 


1 if i = j 
0 if i / j 


(65) 


9A i_ 

ax A 

do 

a o 

3m _ 

3m 

ax. 

ax j 

9A i 

ax. 

axj 

ax i 

i 

1 


( 66 ) 
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TABLE I. - PROBLEM CONSTANTS 


Problem variable 

Assumed value 

3 2 

Gravitational constant of Sun, km /sec 

1. 32715445X10 11 

Gravitational constant of Earth, km^/sec^ 

3. 986032X10 5 

Astonomical units (AU), m 

1. 49599X10 11 

Earth radius, m 

6 378 165 

Initial altitude above surface of Earth, km 

185 

Period of Earth about Sun, days 

365.256 

Final radius, r, m 

2. 278X10 11 

Final velocity, V, m/sec 

24 100 

Specific impulse, sec 

5000 


TABLE II. - SAMPLE CASES FOR TRANSFERS FROM 185-KILOMETER 
ALTITUDE PERIGEE AT EARTH TO ORBIT OF MARS 


Angle between Earth at time zero and final rocket position, 225°. 


Case 

Total trip 
time, 
days 

Type of 
initial orbit 
Earth 

Characteristic 

velocity 

increment, 

AV, 

m/sec 

Initial thrust- 
weight ratio, 
F/W 

Propellant 

fraction 

1 

275 

Parabolic 

8 074 

1 

o 

t-H 

0. 15185 

2 

275 

Circular 

15 853 

10' 4 

. 27626 

3 

275 

Parabolic 

11 318 

0. 56xl0 -4 

. 20613 

4 

240 

Parabolic 

20 390 

10 -4 

. 34022 





Figure 1. - Coordinate system definition. 



X, Earth radii 


Figure 2. - Earth escape path for three-body 
solution. Initial Earth orbit, parabolic 
(case 1); initial thrust-weight ratio, 10 . 


-Thrust vector 



(a) Initial Earth orbit, parabolic (case 1). 


End of tangential thrust; 
start of optimal thrust-x 20 


60 -40 

X, Earth radii 


Figure 3. - Earth escape path for three- body solution (last two turns). Initial Earth 

„4 

orbit, circular (case 2); initial thrust-weight ratio. 10 . 


Y, AU 1. 5 1 — 


’-Earth escape point 



-1.5 j -1.0 -.5 0 .5 1.0 


Coast 

Thrust vector 


(b) Initial Earth orbit, circular (case 2). 


Figure 4. - Heliocentric trajectories for three-body solutions. Initial 
thrust-weight ratio, 10~ 4 . 


Angle of attack, deg 


CO 

CO 



(a) Initial Earth orbit, parabolic (case 1). 



(b) Initial Earth orbit, circular (case 2). 

Figure 5. - Thrust angle program for three-body solution. Initial thrust-weight 
ratio, 10 . 



Earth radii 

Figure 6. - Effect of using tangential thrust instead of 
optimal thrust at beginning of flight for three-body 
solution. Initial thrust-weight ratio, 10"**. 


NHSMm*-' 



Ratio of two- body propellant to 



Figure 7. - Propellant comparison between three-body solution and 
two-body solution with tangential thrust in planetocentric space. 



Patch radius, r g , Earth radii 


Figure 8. - Propellant comparison between three-body 
solution and two-body solution with coasting in planeto- 
centric space. Initial Earth orbit, parabolic. 
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